Cooper pair formation in trapped atomic Fermi gases 
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We apply the closed time-path formalism to evaluate the dynamics of the BCS transition to the 
superfluid state in trapped atomic 6 Li. We find that the Fokker-Planck equation for the probability 
distribution of the order parameter is, sufficiently close to the critical temperature, identical to the 
equation that describes the switching on of a single-mode laser. 
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After the successful experiments to trap and cool 
atomic gases of the bosonic alkalies 87 Rb, 7 Li and 23 Na 
below the critical temperature for Bose-Einstein 
condensation, the field of atomic physics has attracted 
a lot of interest. Besides the ongoing experiments on the 
properties of the Bose condensate in these gases, several 
groups have started experiments in order to observe also 
quantum degeneracy effects in fermionic gases. A partic- 
ularly promising candidate in this respect is the fermionic 
isotope 6 Li which is predicted to make a BCS transition 
to a superfluid state. Indeed, the triplet s-wave scattering 
length a s of this atom is enormously large and negative 
Q , resulting in a critical temperature of the same order 
of those of the BEC experiments ]5| . A requirement for 
this relatively high critical temperature is that a mixture 
of two different hyperfine levels of the 6 Li atom must 
be trapped. Depending on the precise hyperfine states 
trapped, this can be done in a magnetic or optical trap. 
In either case, a study of the dynamics of this BCS phase 
transition offers the exciting opportunity to observe the 
evolution of a spontaneous breaking of symmetry under 
almost ideal conditions. 

In case of a magnetically trapped gas mixture, a draw- 
back of the large s-wave scattering length is the corre- 
sponding large decay rate and resulting limited lifetime 
of the gas . Fortunately, this decay can be suppressed 
by applying a magnetic bias field. For example, the life- 
time of a gas with a typical density of 2 x 10 12 cm~ 3 
atoms per hyperfine state, is of the order of 1 second in 
a magnetic bias field of about 7 T, whereas it is of the 
order of only 0.2 ms at 0.1 T. Evidently, the lifetime of 
the gas must be longer than the time scale on which the 
formation of Cooper pairs occurs. In previous work || , 
we estimated on basis of mean-field theory the nucleation 
time to be 0(h/kBT c ). For a homogeneous mixture with 
a density of 2 x 10 12 cm~ 3 in each hyperfine level, the 
critical temperature was calculated to be 37 nK, resulting 
in a nucleation time of 0.2 ms. Since low magnetic fields 
are experimentally much more convenient, it is also from 
this point of view sensible to make a precise calculation 
of the nucleation time of the BCS transition. The aim of 
the present paper is therefore to evaluate the dynamics 
of the BCS phase transition, from which, as a byproduct, 
the time scale for the formation of Cooper pairs can be 
extracted. 



The dynamics of the phase transition can be de- 
scribed by using the closed time-path formalism [Q. 
This can easily be understood from an analogy in quan- 
tum mechanics where the probability amplitude for a 
particle moving from a position xo at time to to an- 
other position x at time t is given by (x, t|xo,io) — 
J 2?[x] exp (iA\x\/K) with _4[x] = J ( * dt'£[x, x] the action 
and C the lagrangian of the particle. The probability dis- 
tribution for being at position x at time t, i.e., P[x, t] = 

J dx (x,t|xo,to)(xo,io|x,£)P[x ,io], is now & iven h ^ a 
path integral over the periodic paths [x(i')] where t' runs 

from t to to (antichronological, lower branch) and back 
from to to t (chronological, upper branch), i.e., over a 
closed time or Keldysh contour . 

In our case, we want to study the dynamics of the 
BCS phase transition. Analogous to the simple quantum 
mechanical description for the position probability of a 
particle, we thus need to derive the time dependence of 
the probability distribution for the order parameter A 
which describes the superfluid phase The start- 

ing point is the probability P[ip, ip*,t] for a gas of atoms 
in two hyperfine states denoted by \a) = | T}> | I) 



P[W,t]= / P[Vp[^]exp(^[^,V*; 



where the functional integral is over the fermionic fields 
tjj a (x,t) and ^/>*(x, f) that are defined on the Keldysh 
contour C with to — > — oo. The action A[ip,tfi*] is given 
by 



A[*l>,ip* 



n 2 v 2 
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Vo 



where Vo < represents a local attractive interatomic 
interaction. Introducing by means of a Hubbard- 
Stratonovich transformation the order parameter A(x, t), 
the action becomes quadratic in the fermionic fields and 
ij) and ip* can easily be integrated out jnj . We thus ar- 
rive at an effective field theory for this order parameter 
A(x, t) that equals 



A[A, A*] = -iftTrUnG" 1 ] + j dt J dx 



jA(x,t)| 2 
Vo 



, (1) 



1 



where the trace runs over the spin, time and spatial 
variables, and the Green's function matrix is defined 
as G(x,t;x',f) = -z(T c [^(x, t)^*(x', t')]). The spinor 

ijj*(x,t) is given by ip*(x,t) = (^p^(x,t) , -0t( x , *)) an d 
Tc is the time ordening operator on the Keldysh contour 
C, taking into account the anti-commuting nature of the 
fields. 

Around the critical temperature T c , the order param- 
eter A is small compared to fc^T, and Eq. ([j]) can be 
expanded in powers of A. To include to lowest order all 
relevant physics, it suffices to take only the terms propor- 
tional to |A| 2 and |A| 4 into account. To calculate these 
quadratic and quartic parts of the effective action for 
A explicitly, we must remember that the time integrals 
within these expressions are over the Keldysh contour, 
and must be rewritten as time integrals over the real axis 
— oo < t' < t only. The Green's function Go(x, t; x', t') 
defined on the contour, must then be decomposed into 
retarded, advanced and Keldysh components. Further- 
more, the order parameter A(x, t) on the contour can 
be decomposed in a classical part <j) and fluctuations £ 
defined on the real time axis as 



A(x,t±) 



(x,i)±-£(x,t) 



where the upper and lower sign refer to the upper and 
lower branch of the Keldysh contour, respectively. The 
important point is that the fluctuations allow the clas- 
sical part of the order parameter to evolve from zero to 
a nonzero value when the system is driven through the 
phase transition, and we can thus indeed describe the 
dynamics of the classical field </> around the critical tem- 
perature by means of this closed time-path method. In 
contrast, this cannot be achieved within the framework 
of mean-field theory, where the order parameter always 
remains zero if it is zero initially. 

Applying subsequently a gradient expansion which as- 
sumes that spatial and temporal variations in the order 
parameter are small on the size of the Cooper pairs, the 
quadratic plus quartic part of the effective action Eq. (|l|) 
up to quadratic order in the fluctuations, becomes 

A[<j>,<i>*,t,e] = J dx J dt'{0(x,iO£ (fl W)f cm') 

+ £(x,t')£ (A) (x, t')4>* (x,0 

+ ^(x,t')i (K) (x,i')r(x,i')}, (2) 

where to lowest non- vanishing order in V and in d/dt the 
retarded and advanced factors are [p|JfO|] 



(x,f) = a + 7V 2 - /?|0(x, t)| 2 ± <A| (3 ) 



k B T dt 



and the Keldysh part is given by 

L( J O(x,t)=^V |. 



Here we considered only the optimal case fi-\ = /ij and 
used particle-hole symmetry around the Fermi level of 
the system. The constants a, (3 and 7 are given by 



a = N I I - 

= N 7 



T 
T 
C(3) 



8 (TTk B T) 2 

7C(3) (hv F \ 
N °8~6^{k^T 



(5) 
(6) 
(7) 



(4) 



Furthermore, N = mkp / (2-K 2 h 2 ) denotes the density of 
states at the Fermi level, the critical temperature obeys 
k B T c ~ 5/3e F exp(-f/A) j§ with A = 2k F \a s \/iT, and 
Up- is the Fermi velocity of the system. 

The last term in the right-hand side of Eq. (0) describes 
the quadratic fluctuations. If one would neglect them, 
and integrate out the (in that case linear) fluctuations, 
one would obtain a (dissipative) non-linear Schrodinger 
equation for the classical part <j> of the order parame- 
ter which would also arise from equilibrium theory. Of 
course, for our non-equilibrium purposes, the quadratic 
fluctuations are important. The Fourier transforms of 
Eqs. (|) and (@), U R ' A \\a, w) and L^(k, w) in principle 
are related through the fluctuation-dissipation theorem 

i W (k,w) = ±2i(2N B (hu) + l)lm{L {K) ' {A) (\<L,uj)] (8) 

where N B (hu) denotes the Bose distribution evaluated 
at Huj. At low frequencies it thus gives L^ K >(\c, u) = 
±2i(2k B T/huj)lm[L < - R ^^(k,u;)}, in agreement with the 
above results. Notice that if we expand the quadratic 
part of the action A in Eq. ([!]) up to lowest non- 
vanishing order in V and d/dt, we only have to ex- 
pand the quartic part of A to zeroth order, i.e. we only 
need to consider the local contributions proportional to 
<^>(x, i)£*(x, t)\(j)(x, t)\ 2 and it's complex conjugate. Since 
the factor f3 has no imaginary part, we conclude that it 
is indeed consistent to neglect Keldysh-like terms that 
would in principle also arise from the quartic part of the 
action A. 

In order to arrive at an equation of motion for the prob- 
ability distribution of the classical value 0, we first inte- 
grate out the fluctuations £ from the action. Before doing 
so, we however need to realize that in real experiments, 
the gas will be trapped in an external potential V ex ^(x) 
generally in the shape of a 3D isotropic harmonic oscilla- 
tor V(x) = l/2mcjQX 2 . As a result, the coefficients a, (3 
and 7 in Eqs. (^) to (f7j) now depend on x, since the Fermi 
energy /i a (x) = fi a — V(x) — (47ra s ?i 2 /ra)n_ Q (x) and 
therefore the densities n a (x) and the critical temperature 
T c (x) are position dependent. Hence, we can make use 
of the following expansion in the effective action Eq. (0) 
<Kx,*) = E„M*¥n(x) and £(x,t) = £„ & nW0n(xj, 
where the normalized wave functions <p n (x) are solutions 
to Schrodingcr's equation [Q 



2 



- 7 (x)V 2 ](/> n (x) = 7V o (x) £ „0 n (x). 



(9) 



It is well-known that close to T c , the spatial region where 
the order parameter becomes non-zero, is small compared 
to the region over which the density changes ||[l2|]. So, 
we can take N (x) ~ Nq(0), «f(x) ~ vp(0), but we must 
expand T c (x) around x = 0, since it depends exponen- 
tially on the coupling constant A(x). 

Considering again only the optimal case that (i[ = /ij 
and using that the density in each hyperfine level around 
x = can be approximated with n Q (x) = 1/2 n t (0)(l — 
(x/i?) 2 ), where n t (0) is the total gas density in the center 
of the trap, one finds from the zero temperature expres- 
sion for the respective chemical potentials that 



R 2 = l 4 [ 47m 4 (0)a + 2 



V3 



2/3 



[«t(o; 



12/3 



where / = \jTi/muiQ. From this expression one can im- 
mediately see that the effect of the mean-field interaction 
in the gas is to contract the gas cloud to the center of the 
trap. Expanding T c (x) around x = 0, one arrives at [1X2] 



a(x) = JVq(O) In 



T c (0) 6 F (x) \ _ J_ / JV (0) 
T e F (0)J A(0) \N (x) 



- 1 



and Eq. (g) turns out to be just the Schrodinger equation 
for a 3D harmonic oscillator with effective frequency 



OJA = 47Ti 



' 2A+1 k B T 
7AC(3) mv F (0)R 



The eigenvalues are given by 
2m 7 (0) ( 3\ 



In 



T c (0) 
T 



(10) 



Notice that depending on the temperature, the eigen- 
values e„ become unstable, i.e., with decreasing tempera- 
ture they can become negative one by one. The tempera- 
ture where £o becomes zero, is the critical temperature T c 
of the inhomogeneous system. It is clear from Eq. jl0| ) 
that due to finite size effects, this critical temperature 
T c is lower than the homogeneous critical temperature 
T c (0) corresponding to the density in the center of the 
trap [l2]| . Now it is important to realize, that there is 
always a region sufficiently close to the critical tempera- 
ture T c of the system, where it is a good approximation 
to take into account only one single mode 



)(x) = 



/ 1 



3/2 



exp 



211 



(11) 



in the expansion for 4> and £ because the higher modes 
4> n remain stable and their influence thus remains small. 
Furthermore, as will be shown explicitly by a numeri- 
cal example later on, the coupling between modes in the 
nonlinear terms of Eqs. (0) and (0) can be neglected 



sufficiently close to the critical temperature. Plugging 
everything into the action A in Eq. (g) , we can first per- 
form the integral over x, and subsequently integrate out 
the fluctuations that are now solely represented by the 
fields bo(t) and bffi). From the resulting action the effec- 
tive hamiltonian for the classical fields ao(t) and aj(t) can 
be obtained by introducing their conjugated momenta in 
the standard way. We find 

HA = ^l2k B Tj^ 7+ ° [ N (0)e 



7riVo(0) 



da da^ da 
d 



+(3'\a Q \ 2 }a Q + —[N (0)e + /?'|a |>o j , 

where 0' = (3 J dx |0 o (x)| 4 = (3/(2tt1 2 a ) 3 / 2 . 

The equation of motion, or Fokker-Planck equation, 
for the probability distribution P is now given by 
Schrodinger's equation 



ifi^-P = H A P 
at 



(12) 



which has an equilibrium solution 



PiflQ, Oq) oc exp 



where the potential V(ao,ciQ) = A r o(0)eo|ao| 2 + 
1/2 /?' | do | 4 - Since (3' > 0, this potential has the shape 
of a parabolic well for T > T c where £o > 0, whereas it 
becomes mexican-hat shaped below the critical temper- 
ature where £ < 0. 

The dynamics of the phase transition can most easily 
be extracted if we introduce dimensionless variables /, 8, 
t, and the constant a according to 



a 



exp i6 



(13) 



2/3' 

ttNq(0)?i 
^/32(k B TW 

k B T 0' 

The Fokker-Planck equation Eq. (|l2] ) can then be rewrit- 
ten in the form 

d f d 2 1 d 2 d 1 

which exactly is the Fokker-Planck equation associated 
with the single-mode laser 13 1. The constant a repre- 
sents the 'pump parameter', and / is equivalent to the 
'laser intensity'. 

The second term on the right-hand side of Eq. ( |l4| ) 
gives rise to phase-diffusion uM, but when the system is 
initially in a state which does not depend on 6, it will re- 
main phase independent, and we can omit this term from 



3 



the Fokker-Planck equation. The stationary solution in 
terms of the dimensionless variables becomes 



P at (I,a) = A/"(a)exp - 



{I-af 



(15) 



where jV(a) is a normalization factor. This is a Gaus- 
sian distribution centered around I — a, where a < for 
T > T c , and a > for T < T c . Note, however, that / 
only assumes positive values, so even above the critical 
temperature, the average (/) is larger than zero. 
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FIG. 1. Evolution of the probability distribution for the 
order parameter in normalized units. The inset shows the 
average value (7), and its variance as a function of time. 

One can model the phase transition by changing the 
constant a suddenly from a negative value oo < to a 
positive value a > at r = 0. Indeed, the process of 
evaporative cooling is very fast, due to the large value 
of the s-wave scattering length and corresponding large 
cross section for elastic collisions. As a consequence, the 
'translationaP degrees of freedom of the atoms in the gas 
are always in thermal equilibrium, but the order parame- 
ter is not and needs more time to adapt. The probability 
distribution P(I, t, a) will therefore evolve from the sta- 
tionary solution P s t(I, ao), see Eq. (|l5|), at r = to the 
stationary solution P s i(I, a) at sufficiently large r. We 
have solved Eq. (|l^) numerically for ao = — 2 and a = +4 
using the methods developed in Ref. p3| . The evolution 
of the probability distribution P(I, r) is plotted for sev- 
eral times in Fig. [|. It is clear that after a dimensionless 
time of 0(1), the order parameter of the phase transition 
has reached its stationary value. This is shown explicitly 
in the inset of Fig. [ll where the average (I)(t) and the 
variance ((I- (I)) 2 )(r) are plotted. 

In order to see to what extent the approximations we 
have made are valid, we finish the paper with an ex- 
plicit numerical example. For a gas mixture with density 
n a (0) = 2 x 10 12 cm -3 , or equivalently n t (0) = 4 x 10 12 
cm -3 , one has for the homogeneous critical temperature 



T c (0) ~ 37 nK ||. If this gas is trapped in a 3D isotropic 
harmonic oscillator with hujo/ks = 7.1 nK the extend 
of the cloud is R ~ 1 1Z at T — 30 nK, whereas Ia — 2Z is 
then indeed much smaller. The critical temperature T c 
according 'to Eq. @ is T c ~ 31 nK. The second mode 
of Eq. ( |l0|) becomes unstable only at T ~ 26 nK. At this 
temperature the pump parameter a ~ 13, which gives a 
first requirement for our theory. However, a more strin- 
gent one is that it is only a good approximation to take 
one single mode into account if we can also neglect the 
nonlinear coupling in the hamiltonian. We must then re- 
quire that the interaction /3|</>(0)| 2 is small compared to 
N (0)huj A . This is true for T > 29.3 nK, from which it 
follows that the approximations are valid if a < 4.3. For 
T = 29.3 nK, the order parameter A(0) /k B T ~ 1. 

In Fig. [I] we therefore have chosen a = 4, and from the 
inset of this figure we conclude that for a total density 
rtt(0) — 4 x 10 12 cm~ 3 , one unit of dimensionless time 
r corresponds to 8.9 ms, which is much larger than the 
fr/fcgT c (0) = 0.2 ms as estimated in our previous work 
||, because Eq. ( |l3| ) clearly shows that the time scale 
depends also on several other relevant parameters. In 
addition, this is caused by a critical slowing down, which 
can be reduced by increasing the final value of the pump 
parameter a, but requires inclusion of more modes in the 
theory. In a future publication we will come back to this 
multi-mode regime. However, within the present theory, 
for a lifetime of the gas of the order of the nucleation 
time, one needs a magnetic bias field of about 0.7 T. A 
better option is to trap and cool the lowest two hyperfine 
states in an optical trap, which has recently been shown 
to be possible for bosons jlj|. In any case, we believe that 
this paper, in combination with Ref. Q, shows that there 
are no theoretical limitations to achieve a BCS transition 
in atomic 6 Li. 
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